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Time-dependent local Green's operator 
and its applications to manganites 
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An algorithm is presented to calculate the electronic local time-dependent Green's operator for 
manganites-related hamiltonians. This algorithm is proved to scale with the number of states 
A'^ in the Hilbert-space to the 1.55 power, is able of parallel implementation, and outperforms 
computationally the Exact Diagonalization (ED) method for clusters larger than 64 sites (using 
parallelization) . This method together with the Monte Carlo (MC) technique is used to derive 
new results for the manganites phase diagram for the spatial dimension D=3 and half-filling on a 
12x12x12 cluster (3456 orbitals). We obtain as a function of an insulating parameter, the sequence 
of ground states given by: ferromagnetic (FM), antiferromagnetic AF-type A, AF-type CE, dimer 
and AF-type G, which are in remarkable agreement with experimental results. 
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I. INTRODUCTION 

The doped manganites are interesting not only because 
of the potential applications of their 'colossal' proper- 
ties, but also because they are fascinating systems to 
study due to the delicate balance of interactions between 
charge, spin, and orbital degrees of freedom. [J 

One of the successful theoretical approaches to the 
study of the manganites phase diagrams, is based 
on effective models that consider the competition be- 
tween double-exchange (DE), superexchange (SE), and 
the electron-phonon interactions 2]. These models and 
techniques have predicted interesting phenomena, like 
nanoscale phase-separationj3^ , influence of disorder in 
metal-insulator transitions 3, and the existence of non- 
trivial magnetic phases Isj, among other phenomena. The 
ground states of such systems and the finite-temperature 
properties were obtained using simulations on clusters 
of spins using the Monte Carlo (MC) technique together 
with the calculation of the electronic energies by means of 
the exact diagonalization (ED) techniquel^ at each sin- 
gle MC step. In order to extract useful information from 
the simulations, one needs to analyze different sizes for 
the clusters, where the CPU (Central Processing Unit) 
memory and time scale as the number of states N in 
the Hilbert space to the third- and fourth-power, respec- 
tively. In 1999, the Truncated Polynomial Expansion 
(TPEM)[^ was devised in order to reduce the CPU time 
scaling with N, but unfortunately at a cost of a large 
prefactor and a detailed comparison with exact case (ED) 
was not shown. 

Moreover, it is of fundamental importance to consider 
the Mn-Cg orbitals dx^_y2 and d^^2_^2 at each siteQ, 
which doubles the size of the matrices respect of the case 
of a single band. The ED operations have to be per- 
formed several thousands times in a typical simulation, 
practically limiting the spatial dimensions of the clusters 
considered to D=l and 2. The reduction in the dimen- 
sionality of the system artificially breaks the degeneracy 
of the orbitals d^2_y2 and d^^2_j.2, modifying the elec- 
tronic properties of the system. 



For all these reasons it is important to find an alter- 
native way to calculate the electronic energy that avoids 
the natural limitations of the ED scheme, allowing to 
expand the size of the clusters and the spatial dimen- 
sionality. An efficient and fast way to calculate the time- 
dependent Green's operator in an effective model con- 
ceived for manganites will be described in this work. This 
method uses the Chebyshev expansion of the time evolu- 
tion operator^ which is an extremely accurate and fast 
way to obtain the dynamics response of quantum many- 
body systems. For a system consisting of 10 spins S=l/2 
with Heisenberg-like interactions, this expansion is three 
orders of magnitude faster than the ED method, and for 
a variety of considered cases, is one of the most efficient 
time-marching schemes to keep track of the time evolu- 
tion of a quantum system, where CPU memory and time 
basically scale linearly with system size.0 

In this work we will focus our attention to non-trivial 
single-particle Hamiltonians, where the Hilbert space 
grows linearly with the size of the system. More gen- 
erally speaking, the presently described method, which 
will be referred as the Dyn method, provides an alterna- 
tive way to calculate the local time-dependent Green's 
operator for the time-dependent Schrodinger equation 
(TDSE). The Chebyshev expansion of the hermitian, 
time- independent, differential Hamiltonian operator H 
will be applied to a model related to manganites, which is 
more efficient than the ED technique for electronic clus- 
ters larger than 64 orbitals (using full parallelization). 
The feasibility of the method will be tested proposing 
fixed spins configuration and comparing the densities of 
states with well known analytic results. Then the ground 
states, electronic properties and phase diagram for D=3 
will be determined by means of the Dyn method together 
with the MC technique and a comparison with experi- 
mental results will be made. 

The text is organized as follows: in Section II the 
Hamiltonian model is described, the formalism regarding 
the calculation of the time-dependent local Green's oper- 
ator is given in Section HI and in Section IV a benchmark 
test is used to compare the CPU time as a function of the 
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lattice sites needed by the ED, TPEM and Dyn methods. 
In Sections V and VI, the cases of one orbital per site are 
considered. For fixed configurations of the spins, the elec- 
tronic density of states are calculated and compared with 
analytical results, and the relative error is discussed. In 
Section VII the Dyn method is used together with the 
Monte Carlo technique to obtain the ground state and 
reproduce the phase diagram in two spatial dimensions 
(D=2), and the results are compared with the MC+ED 
technique. In Section VIII the MG+Dyn technique is 
used to obtain the phase diagram for filling x=l/2 and 
D=3 on a 12x12x12 cluster, and the new results are dis- 
cussed. Finally, Section IX contains a discussion of the 
the main results and general conclusions. 



II. MANGANITE EFFECTIVE HAMILTONIAN 

We will consider throughout this work the cases of (A) 
one and (B) two orbitals per site, together with the spa- 
tial dimensions, D=l, 2 and 3. The Hamiltonian model 
to be consideredj5j is quadratic in fermionic operators, 

ia-yy'cr i 

+ -'^AF ^ Si • Sj + A ^{QiiPi + Q2iTxi + QsiTzi) 

(iJ) i 
+ il/2)J2if3rQli + Qli + Ql), (1) 



where the degeneracy is recovered when the hopping pro- 
cesses along the a=z direction are included, that is for 
tl.=tlb=tL^O and tl^^iU.t^/S. 

In the second term of Eq. (1), the Hund constant 
Jh(>0) couples the spin 8;= 4-yui'^'^i'^2d'i'f,y2 
(cr=Pauli matrices) of the Cg electrons with the localized 
t2g-spin Si. The constant Jh is here considered as infinite 
or very large, and the model is drastically simplified. Q 

The third term is the AF coupling Jaf between NN t2g 
spins. The fourth term couples eg electrons and MnOe 
octahedral distortions, A is a dimensionless coupling con- 
stant, Qii is the breathing-mode distortion, and Q2i and 
Qsi are, respectively, the {x^—y"^)- and (3z^— r^)-type 
Jahn- Teller (JT) mode distortions, p-^— J2-f a '^l-yrr'^n'r^ 

^xi= J2ai4i,adiha +4^^^di^a), and r^i= Ea(4a^^ia<T 

— djjj^dibcr)- The fifth term is the quadratic potential for 
adiabatic distortions and is the ratio of spring con- 
stants for breathing- and JT-modes, which for mangan- 
ites is approximately given by «2.[8| 

In the (A) case, one spherical s orbital per site will 
be considered and the orbital indexes in Eq. (1) can 
be dropped, that is 7 = 7'. The hoppings amplitudes 
for the different spatial cases are defined for (A) D=l as 
t^^tijth] for (A) D=2 as V^^t'y=Ujth and for (A) D=3 
as t^^ty=t^=tijtfi- For the case (A) the electron-phonon 
interaction and the oxygen degrees of freedom are not 
considered, that is (5ii=Q2i=Q3i=A=/3^=0. This case 
will allow to compare the density of states obtained with 
the Dyn method and well-known analytical results. (2j 



For the case (B), the operators diao- (rfibo-) represent 
the annihilation of an Cg-electron with spin cr, in the 
dx2-y2 (d^z^-r^) orbital at site i, and ^ is the vector 
connecting nearest-neighboring (NN) sites. The first 
term in H is the NN hopping of Cg electrons with 
amplitude t^^, between 7- and 7'-orbitals. For the 
(B) D=2 case, the hopping amplitude along the |- 
direction is given by: t^^=-^/3t^^^^-V^t^^=5tl^=Ujth, 
and t^^—^/3t^^=^/3t^^=3i^^=tijth- The parameter th is 
the hopping transfer integral and will from now on the 
energy unit of the model, and the parameter tij is a com- 
plex scalar that depends on the relative orientation of 
neighboring localized spins S^ and Sj, assumed classical 
with |S|=1, and characterized by the polar and azimutal 
angles, 6i and (/jjQ 

Uj = cos(6ii/2) cos(6'j/2) + e-^('^.-'^^) sin(6l,/2) sm{0j/2) 

. (2) 

when the approximation of very large Jh is used. It 
should be kept in mind that for A = the original de- 
generacy of the dx2-y2 and d^z'^-r^ orbitals is explicitly 
broken in the D=2 case. The hopping process along the 
X- and y- directions favor the occupation of the d^2_y2 
over the d^z^_j.2 orbitals in an important range of fillings. 
For that reason it is important to consider the D=3 case. 



III. TIME-DEPENDENT LOCAL GREEN'S 
OPERATOR 

Without loss of generality, the bra and ket notation will 
be used from now on. Although the formalism given in 
this work is oriented to analyze a model related to man- 
ganites (Eq.l), the results are valid for any hermitian, 
time- independent, differential operator H expressed in 
matrix form, which possesses a complete set of eigenfunc- 
tions {| and eigenvalues {7„}, satisfying equa- 

tions of the form, 

H\K)^ 7„ I 0„) (3) 

In case the set of eigenfunctions and eigenvalues 
are known, the Green's operator for a one-particle 
HamiltonianQ G{t) can be calculated as, 

N 

G(t) = -le-^""^ = ^ e-^"*/'-' I 0„) (0„ I (4) 
11=1 

where the Planck constant will be set to ft = 1. We 
can define G+(i) = 9 {t) G{t), where 6 (t) is the Heavi- 
side step function. After Fourier-transforming G~^{t), the 
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frecuency-dependent Green's operator is obtained, 



where 



+ 00 



''^G+{t)dt 



27r 



(5) 



which in turn, allows to obtain the electronic density of 
states, 



p{uj) 



-Im [Tr {G+{uj))] 



(6) 



However, in order to evaluate (Eq.1-6) the knowledge 
of the eigenfunctions is needed, which is a difficult task 
even for non-trivial single particle Hamiltonians. In this 
work a novel way is proposed to calculate this local time- 
dependent Green's operator, without the explicit knowl- 
edge of the set of eigenfunctions {| </>„)}. It will be shown 
later (Fig. 1) that this algorithm outperforms the ED 
technique for lattices larger than 64 sites (using paral- 
lelization), without providing the information regarding 
the eigenfunctions {| The local time-dependent 

Green's operator must be expressed in terms of the lo- 
cal basis states {| ?^)} instead of the eigenfunction's set 
{I 4'n)} as 



G{t) 



N 



iHt 



(7) 



and then follow the steps of Eqs. 5-7 to obtain p{lo). 
It has been shown by Dobrovitski et al.Q that using a 
Chebyshev expansion method can be an extremely pre- 
cise and fast way to evaluate the time evolution opera- 
tor. The time-dependent Green's functions are matrix 
elements of the operator G, and are evaluated as 



Jkir) 



e-"*Tk{x)dx 



(10) 



are the fc-order Bessel function of the first kind and 
Tk{x) are the /c-order Chebyshev polynomials of the first 
kind, given by: Tk{x) = arccos(A; cos(a;)). The vectors 
I i^k) are calculated following the Ghebyshev recursion 
expression: | vq) ~ 1- | ig), | v\) — X- \ jq), and 
\ v,,) = X \ Vk-i)- I Vk-2) (for k ^ 2). Since the 
value of a Bessel function decreases as Jk (t) « (jlk)^ , 
the truncation of the series to the order K^ax leads to an 
error that decreases exponentially with Kmax- In prac- 
tice, holding terms of the order Kmax ~ 1.5t is enough 
to get an accuracy of 10^^ in the wave function. 



IV. BENCHMARK TESTING 

The performance of the Dyn algorithm was bench- 
mark tested on AMD Athlon 2500-h (1.85Ghz) proces- 
sors against the ED algorithm and TPEmT^] on square 
random hermitian matrices. In Fig. 1 are shown the wall- 
clock computational times t required to compute Eq. (5) 
in one Monte Carlo step per site as a function of the rank 
of the matrix H, Rh, using the ZHEEV - LAPAGK 
library (ED algorithm, triangles), TPEM (squares) and 
the Dyn algorithm (circles). The computational times 
t were fitted in the range 100 < Rh < 1000, obtain- 
ing the following dependencies: logitEo) = —5.47 + 
3.511og(i?ff), \og{tTPEM) = -4.59 -h 2.05 log(i?ff) and 
log(iDyn) = -3.15 -I- 1.551og(i?ff). For the TPEM and 
Dyn cases the full parallelization capability was used and 
Kmax = 60 was considered. 



N 



10000- 



G{io,jo;t) = -i^Uo \ n){n 



iHt 



I *o> (8) 



where | zq) and | jo) G {I and the diagonal ele- 
ments G{io,io;t) are relevant for determining p{lu). The 
quantity | G{io,jo;t) ^ is the probability of creating the 
particle at site io and detecting it at site jo, at time a t 
later, and a decaying behavior with time is expected. 

In order to carry out this expansion, it is first nec- 
essary to normalize iJ, by a value | 7,jjax I equal or 
higher than the highest eigenvalue in absolute value: 
X = H/ I 7„,ax I- the Ht phase must be kept con- 
stant, the time t is normalized too, hy t = t\ 7max I-E3 

Now the expansion of the normalized Hamiltonian X 
operator at time r is, 



(jo I e-^"' I zo) 



Jo(t) I vo)+2Y,Jk{T) 

fc=i 



I'k) 



(9) 
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FIG. 1: Computational time t required to evaluate Eq. (5) 
in one Monte Carlo step per site as a function of the rank 
of the matrix H, Rh, using ED {ZHEEV - LAPACK li- 
brary, triangles), TPEM (squares) and Dyn (circles). Full 
parallelization capability was used for the TPEM and Dyn 
codes. 
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We can observe that for Rh ^ 64, the present method 
and TPEM outperform the ED technique. For matri- 
ces smaller than 144x144, the TPEM technique is faster 
than Dyn, but this situation is reverted for matrices of 
size Rh > 289. This extraordinary performance can 
be understood intuitively by the following reasons: 1) 
The efficience of the recursive relation of the | i^k) vec- 
tors, 2) The functions JkW) are relatively inexpensive 
to obtain computationally [l3|, 3) No explicit multiplica- 
tion between the T]^{x) operators and the vectors | zq) is 
necessary. 4) The use of the parallelization capability. 
If only one CPU is used (no parallelization) in the Dyn 
(TPEM) algorithm, an increase in a factor of Rh is ob- 
tained in the time consumption toyn {trPEAi)- In this 
case Dyn outperforms ED only for Rh ^ 900. 



CASE A: ONE ORBITAL PER SITE AND FM 
CONFIGURATION 
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FIG. 2; The probability | G{io, jo',T) p as a function of the 
site jo, when a particle is created at time r = 0, in the middle 
of a 1000-site chain (PBC), and measmed at a time r — 1000 
later (a), and at a time r — 3000 (b). The arrows show the 
direction of movement of the wave- front. 



The Dyn method can be directly tested with the an- 
alytical results for the density of states in cubic lattices 
for D=l, 2 and 3. This can be done by considering one 
orbital per site for Eq. (1), and the case where di and 
tp.^ are constant for all i, that is the ferromagnetic (FM) 
phase where all the spins Si are parallel. The effect of 
the phonons will not be considered and H is normalized 
by the value | j^^^ |= 2Dth. 

The Dyn method basically is an efficient way to keep 
track of the time evolution of a wave-packet. In Fig. 2 
can be seen the probability | G{iQ,jo;T) p as a function 
of the site jo, when a particle is created at time r = 0, in 
the 500th site of a 1000-site chain and destroyed at a time 
r later, considering periodic boundary conditions (PBC). 
In Fig. 2(a) the snapshot was taken at a time r = 1000, 
where the wave-front is moving away from the zq site, 
indicated by the arrows. In the case 2(b) the wave- front 
has crossed the boundaries of the chain, and is moving 
towards the center of the chain, at a time r = 3600. [T^ 

It is possible to recover results for the infinite-size limit, 
provided the maximum time while the simulation Tgim is 
carried out is less than the time that the 'wave front' r^f 
reaches back to the site io, Tsim < ^wf- In the opposite 
limit, Tsirn ^ T^f, this method obtains the discrete res- 
olution of the density of states spectra corresponding to 
a finite-lattice model. 

Keeping track of the wave amplitude at the site iq, 
will allow us to study the behavior of G(io, io \ t) for the 
different spatial dimensions. In Fig. 3 it is shown the 
time-dependence of /m [G(io, zo; t)] for the D=l (black 
squares, 1000 sites), D=2 (red circles, 100x100 sites) and 
D=3 (green triangles, 50x50x50 sites) cases. The corre- 
sponding real component of G{io,io;T) is zero. 

After Fourier-transforming the value of G(io, zo! t), the 
corresponding densities of states p (w) are obtained and 
shown in Fig. 4, and are in remarkably good agreement 
with analytical densities of states. The simulations are 
carried out for the case Tgim < T^f, where an approxi- 
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FIG. 3: Time-dependence of Im[G{io,io',T)] for the D=l 
(black squares), T)—2 (red circles) and D=3 (green triangles) 
cases, for lattice sites 1000, 100x100, and 50x50x50, respec- 
tively. The simulation is carried out up to the cut-ofF time 
Taim- At time t — 0, /m [G(zo, zo; 0)] =1, for D=l, 2 and 3. 



mation to p (w) for an infinite system is obtained. 

In order to compare the discrete nature of the density 
of states of finite clusters, a comparison of p (w) for an 8x8 
cluster of sites is given in Fig. 5, with Tgim— 800 (black 
line), and Tgim— 8000 (red line). The values picked for 
Tsim correspond to the case Tgim ^ T^f, and the greater 
the value of Tsi„i the higher and narrower are the peaks 
obtained for p{uj), for that reason the distributions are 
not plotted on the same vertical scale. The distribution 
of eigenvalues obtained with ED is as follows (in units 
ofth): ±4, ±3.41421, ±2.82843, ±2, ±1.4142, ±0.58579, 
and 0. The eigenvalues ±4 are not degenerate, and the 
rest of them are four-fold degenerate with the exception 
of the eigenvalues ±1.4142 (8 times), and the eigenvalue 
(14 times). The peak positions obtained with Dyn are in 
good agreement with the ED results, where the precision 
is proportional to Tsim ■ 



FIG. 4: Frequency-dependence of p (u)) for the FM case in 
ID (black squares), 2D (red circles) and 3D (green triangles). 




FIG. 5: Electronic density p (u) of states for a 8x8 cluster 
for the FM case in the energy interval [-4.254^,-1. 75th], ob- 
tained using the Dyn method, with Tsim=9>QQ (black line) 
and r aim =8000 (red line). Details of the complete spectra 
found following the ED method are described in the text. 



The finiteness of the time simulation Tsim causes an 
oscillatory dependence in p (w) , and for values of lo close 
to the eigenvalues, p (w) can have even small negative 
values. This behavior is commonly known as 'overshoot- 
ing' or Gibbs oscillations in the Fourier transformation 
context and these effects are treated here following the 
Kernel Polynomial approximation TH]. The peaks distri- 
butions obtained with Dyn have a finite-frecuency width 
Aoj « (4Z?7r) I Tsim-, where 4Z? is the electronic band- 
width in units of th- 

The total electronic energy is obtained by integration 
of p (iS) , and in Fig. 6 its dependence vs. chemical poten- 
tial p is shown, where the full (red squares) and dashed 
lines (black triangles) corresponds to the electronic en- 
ergy obtained with the ED {Dyn) method at inverse tem- 
peratures /3=10 and 100 respectively, for Tsim=400. For 
the depicted w— range it can be observed at low tempera- 
tures two plateaus centered around the eigenvalues —At,h 
and — 3,41421t/j, obtained after integrating the delta- 



FIG. 6: Electronic Energy vs. chemical potential. The 
full (red squares) and dashed lines (black triangles) corre- 
spond to the electronic energy obtained with the ED (Dyn) 
method at inverse temperatures /3=10 and 100 respectively, 
for rsim=400. 



shaped density of states. 



10 



2 1 

"8 0.1 

03 

> 0.01 





P 




1 


-0- 


10 


--A-- 


100 


V 


1000 



100 1000 

lattice sites 

FIG. 7: Averaged error of electronic energies as a function 
of Tsim, for values of /3=1 (squares), 10 (circles), ICQ (up 
triangles), and 1000 (down triangles). 

We define the relative averaged error of the Dyn 
method in the electronic energy, for p G [— 2Dt/i, +2Dth] 
by, 



averaged error — 



■ E 

1=1 



{EED,i — Eoyr. 



EeD.: 



where the electronic energies Eed,! and Eoyn.i were cal- 
culated using both methods for Nr^^^— 64:00 equispaced 
frequencies in the interval [—2Dth,+2Dtii], and plot- 
ted in Fig. 7 as a function of Tsim, for values of 13—1 
(squares), 10 (circles), 100 (up triangles), and 1000 (down 
triangles). We can observe that the averaged error is 
larger when /3 increases, and decreases when Tsim is in- 
creased, but in a non-monotonic way, and it diverges for 
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values of /i out of the band. In the present work, the Dyn 
method wiU be used together with the MC technique for 
filling x=l/2, where the averaged error is less than 0.1%. 
High precision schemes can be achieved for the case when 
Tsim/Tyjf = i > 1 and L < mod{Tsim,Twf), where L is 
an integer number and mod{) is the modulus operation. 



VI. CASE A: ONE ORBITAL PER SITE AND 
PM CONFIGURATION 

The case of a random distribution for 6i and ipf 
Oi G [0, tt] and (p^ G [0, 2tt] is relevant in the present 
model of Eq. (1) as it represents a particular statistical 
realization of the paramagnetic (PM) phase, and it will 
be used as another test for the Dyn method. In Fig. 8 it 
is shown the time dependence of Im [G{io, Iq; t)] for the 
FM (full black line) and the PM (blue open circles) cases. 
In the last case, a random iq site in a D=2, iV=6400 sites 
lattice (80x80) was chosen. We can see an irregular time 
dependence, with an average characteristic period longer 
than the FM case. 
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FIG. 9: (a) Local density of states (uj) for a random spin 
configuration {(fi^,9i} on a 80x80 lattice, (b) Total density 
of states p{lo) = jr J^^^^i Pigi^) (A^=6400) for the PM (red 
line) and the PM (black line) phases. 



VII. CASE B: TWO ORBITALS PER SITE AND 
THE MC ALGORITHM FOR GROUND STATE 



FIG. 8: Time dependence of Im [G{io,io; t)] for the FM (full 
black line) and the PM (blue open circles) cases, for a random- 
picked io site in the 80x80 square lattice. 

The corresponding PM local electronic density Pj^^ (cj) 
is depicted in Fig. 9(a), for the same random-chosen site 
iQ. The total density of states can be obtained after aver- 
aging the local density of states p{uj) = jj J2fo=i Pioi^)^ 
for the particular random {ip^,9i} configuration consid- 
ered (Fig. 9(b), red hue). The disorder in the hopping 
distribution in the PM phase reduces the FM band- width 
(Fig. 9(b), black line) by a factor 1/^/2. The quantity 
p{uj) can be obtained analytically for the PM case after 
considering this reduction of the FM band-width, but 
the ED results are not shown here, since the CPU time 
needed is about 6 orders magnitude longer that the cor- 
responding using the Dyn method. 

As the calculation of Pig{oj) is completely indepen- 
dent of the site considered, a speed-up of the algo- 
rithm should be achievable by means of a computational 
parallel implementation. |16|| 



The full expression of Eq.(l) is taken into account, to- 
gether with the consideration of two orbitals per site. 
The ground state is sampled through the calculation of 
the partition function, where the electronic component is 
given by: 

/+ao 
p(co)(l + e-^("-^^)dL^ (11) 
-oo 

obtained for each fixed proposed configuration the angles 
and displacements coordinates {Oi, (p^,Ux.^,Uy^i,Uz^i} for 
each site, and /3 is the inverse temperature. The model is 
analyzed primarily using a classical Monte Carlo (MC) 
procedure for the localized spins and phonons, in con- 
junction with the ED and Dyn methods for the elec- 
tronic matrix. This last part of the process corresponds 
to the solution of the single-electron problem with hop- 
pings determined by the localized spin and phonon con- 
figuration. The resulting electronic density is then filled 
with the number of electrons to be studied, that is, the 
simulations are carried out in the canonical ensemble. 
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A quantitative comparison of the electronic energies 
between both methods is given in Fig. 10, where the 
total averaged energies per site obtained using the ED 
(black triangles) and Dyn (red circles) techniques, as a 
function of the parameter Jaf /t- The Monte Carlo simu- 
lations were performed on a 4x4 cluster, at a temperature 
r/t=0.025, A=1.0, (3^=2 and filling x=l/2 (one particle 
every two sites). The vertical dashed lines represents ap- 
proximately the values of JAp/t where the energy level 
crossings occur between the FM, the AF-type CE, and 
the AF-type G phases. The phase diagram for D=2 was 
already obtained using the ED method^, and the Dyn 
method is in excellent quantitative agreement with these 
results. 




FIG. 10: Total energies per site obtained for the effective 
manganites model (Eq. 1) using the ED (black triangles) and 
Dyn (red circles) techniques, as a function of the parameter 
JAp/t. The Monte Carlo simulations were performed on a 
4x4 cluster, at a temperature T/f=0.025, A=1.0, = 2 and 
x=l/2. The dashed lines represents approximately the JAp/t 
where the energy level crossings occur between the FM, the 
AF-type CE, and the AF-type G phases. 



MCS/S. The magnetic character of the different ground 
states were analyzed by means of Spin Structure Fac- 
tor S'(q) =l/A^2X;i j(Si.Sj)e'i ('''-''j). For the 12x12x12 
cluster, the q values considered are: (/7r/6,/c7r/6,m7r/6), 
where k and m are integers values in the interval 
< Z,fc,m < 6. 
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VIII. RESULTS FOR D=3 

Novel results for the phase diagram x=l/2, A=0.5, 
D=3 at T/th^Om were obtained following MC+Dyn 
simulations on a 12x12x12 cluster (3456 orbitals). 
The maximum number of momenta considered was 
Kmax=^^^ or 200, and the number of Monte Carlo steps 
per site (MCS/S) was typically taken as 5000. Full par- 
allelization of the algorithm was performed, where tipi- 
cally 288 cpus where dedicated to compute a single task. 
Updates of the spin and phononic {9i,ip^,Ux.i,Uy^i,Uz^i\ 
configurations were accepted or rejected according to 
the Metropolis algorithm. The simulations start in 
most of the cases with random initial configurations, 
but for the A and CE phases the convergence is very 
slow, specially close to the energy crossovers. A speed 
up of the convergence was realized by fixing the cor- 
responding spin configurations, and testing the sta- 
bility of the proposed ground state as a function of 



FIG. 11: (a) S(q) vs. MCS/S for JAF/th=0.3, x=l/2, 
A=0.5, D=3 at r/tfe=0.02 on a 12x12x12 cluster and q = 
(tt, 7r,7r/2) (black squares), q = (7r,7r/2,7r) (red circles) and 
q = (7r/2,7r,7r) (green triangles). The ground state corre- 
sponding to this setup of parameters is the dimer phase, (b) 
Chemical potential (black squares), electronic (red circles), 
superexchange (green up triangles), elastic (blue down trian- 
gles) and total energies (cyan line) vs. MCS/S. 

In Fig. 11(a) it is shown S'(q) vs. MCS/S for 
Jaf /th=^-i and q = (7r,7r,7r/2) (black squares), q — 
(tt, 7r/2,7r) (red circles) and q = (7r/2,7r, tt) (green trian- 
gles) . The ground state corresponding to this setup of pa- 
rameters is the dimer phase, which consists of pairs spins 
tt: aligned antiferromagnetically between them. This 
phase is evidenced by a single peak 5(q) =1/2 at one of 
these q values chosen. For MCS/S<150 the dimer corre- 
lations between neighboring spins are pointing in all the 
three spatial directions. For MCS/S>200, the system 
breaks the spatial isotropy, and all the cluster consist of 
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dimers aligned in the z-direction, for this particular ex- 
ample, and for MCS/S>400 the system has reached ther- 
mal equilibrium. In Fig. 11(b) it is shown the Chemical 
potential fj, (black squares), electronic (red circles), su- 
perexchange (green up triangles), elastic (blue down tri- 
angles) and total energies (cyan line) vs. MCS/S. The 
Chemical potential is calculated autoconsistently at 
each Monte Carlo step to fix the electronic density to 
X — 0.5 ± 0.00001. The electronic energy is given by the 
first and fourth terms in Eq. (1), the superexchange en- 
ergy involves the second term and the elastic energy is 
given by the fifth term. 

The total energies per site of the ground states are 
plotted as a function of the parameter JAp/th in Fig- 
12, under the same conditions used in Figs. 11. As the 
JAp/th parameter is increased, the AF-insulating char- 
acter of the ground states obtained increases. Compar- 
ing this picture with the phase diagrams for D=2 (Fig. 
10), we notice the appearance of the A and dimer phases 
for the case D=3. The A-phase, which consists of FM 
planes aligned antiferromagnetically between them, and 
the FM-phase are degenerate for D=2, but for the case 
D=3 their corresponding total energies are different, i.e. 
the degeneracy is removed. The dimer phase, which was 
obtained in D=l and 20 considering one orbital per site, 
is not a ground state in T)—2 with two orbitals per site 
due to the explicitly breaking of the symmetry between 
dx2_y2 and d^^2_j,2 orbitals. 

Increasing the insulating parameter Jaf / th the follow- 
ing sequence of phases is obtained: FM, A, CE, dimer 
and G. It is worth to note that, with the exception of 
the dimer phase, these phases have been observed exper- 
imentally in half-doped manganites,|17|, , Ifl. 20, 21] 
following approximately the same sequence. |2^ 

For values of JAp/th in the interval [0,0.125] the 
ground state of the system corresponds to the metal- 
lic FM phase. The orbital configuration of this phase 
is highly degenerate, but analyzing snapshots at fi- 
nite temperature orbital states of the form cos(8) | 
dx2-y2)+sm(Q) I d^z^-r^) with random values of Q in 
the range < 8 < 27r are observed, namely the orbital 
state of this phase is completely disordered with an ho- 
mogeneous distribution of the charges. 

The A-phase is found for values of JAp/th in the in- 
terval [0.125,0.16]. The orbital distribution corresponds 
to a majority occupation of the planar dx2_y2 orbitals, 
which favors a metallic in-plane conductivity. At the 
same time this orbital arrangement favors the AF align- 
ment between neighboring FM planes, where the hop- 
ping process inter-planes is strongly suppressed. This 
phase has been observed in the Pro.sSro.sMnOa^] and 
Ndo.sSro.sMnOa^] compounds, in this last case in co- 
existence with the CE-phase. 

For values of JAp/th in the approximate interval 
[0.16,0.26] the ground state of the system corresponds 
to the CE-phase. T9] This magnetic phase was mea- 
sured experimentally in the systems Lan t^Can t;MnO-i|20j]. 
Ndo.sSro.sMnOaiia], and Fro.sCao.sMnOai^ljj. Even for 
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FIG. 12: Total energies per site of the ground states as a 
function of the parameter Jaf /th for the same conditions as 
in Fig. 11. 



A=0, the character of this phase is insula ting and was 
analyzed in several theoretical worksP.l2Ml24j. were the 
details of the spin, charge and orbital configuration can 
be seen. [23 

The dimer phase, which is the ground state of the 
model in JAp/th G [026, 0.66] for D=3, was also obtained 
in D=l and 20. The existence of FM dimers (Zener 
Polarons) along the zigzag chains of the CE-phase has 
being observed in Pro.6Cao.4Mn03'26'] compounds, to- 
gether with a dimerization of the Mn-Mn distances. More 
complex interactions should be included in the present 
model, like the displacement of the Mn cations together 
with the distance dependence of the couplings Jap and 
th in oder to reproduce this phase. For the pairs IT di- 
rected in the ^ = x, y or z spatial directions, the orbital 
configuration of the dimers phase corresponds to orbitals 

Finally, for values JAp/th ^ 0.66 the insulating G 
phase is obtained. The range of parameters where the 
G and dimer phases are present corresponds to mangan- 
ites with very small band-width and sofar has not been 
experimentally observed up to date. 



IX. CONCLUSIONS 

In this work it was described an efficient high- 
performance algorithm to calculate the local time- 
dependent Green's operator for general hermitian, time- 
independent, linear differential operator H , expressed in 
matrix form. The time-dependent local Green's operator 
is expressed as a finite series of the Chebyshev polyno- 
mials of the normalized H operator, multiplied by the 
corresponding Bessel function of the first kindQ- This 
approach, together with the fact that in the local ba- 
sis no multiplication of matrices by vectors is needed to 
keep track of the time-evolution of an initial state, leads 
to an algorithm where CPU time and memory scale lin- 
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early with the number of states in the basis of the Hilbert 
space. 

The total energies of an effective Hamiltonian for man- 
ganites were compared using the Exact Diagonalization 
technique and the novel Dyn method. The results ob- 
tained with both ED and Dyn methods were in general 
good agreement, and for lattices larger than 289 sites 
the Dyn method outperforms the ED and TPEM tech- 
niques. A parallelization of the algorithm is possible, 
with a speed-up factor close to the number of processors 
involved. The Dyn method could become the method of 
choice to expand the present computational limitations 
set by the ED method, and in particular to study the 
model for D=3. 

It was shown that it is possible to recover results for the 
infinite-size limit (14], provided the maximum time of the 
simulation on a finite-lattice is less than the time that 
the 'wave front' reaches back to the site ia, t < Tsim- 
The discreteness nature of the spectra is obtained in the 
opposite limit, r ^ Tsim- 



The Dyn method used together with the Monte Carlo 
technique was used to obtain new results for the mangan- 
ites phase diagram for D=3 and half- filling in a 12x12x12 
cluster (3456 orbitals) . As a function of an insulating pa- 
rameter we found the following sequence of ground states: 
FM, AF-type A, AF-type CE, dimer and AF-type G. 
These phases are in remarkable agreement with experi- 
mental results in half-doped manganites.pTl FlS. 19^2^ 
l2lj | which follow approximately the same sequence. [2^| 
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